

clear all; close all; dbstop if error;
set(groot,'defaultLineLineWidth',2) ;  


%%% options

    % simulations: saving, loading or re-running
options.Run      = 0 ;
options.Save     = 0 ;
options.Load     = 1 ;

    % whether to include human capital accumulation or not in simulations
options.Learning = 'Learning' ;


%%% custom colors

MyOrange = [ 231,101,45 ] ;
MyOrange = MyOrange / 255 ;
MyBlue =  [ 60,106,209 ] ;
MyBlue = MyBlue / 255 ;
DefGreen = [ 0.4660    0.6740    0.1880] ;
OliveGreen = [53 98 26 ] / 255 ;


%%% solve the model under different scenarios

if options.Run 
    

    %% load parameters estimates
    load('estimates/ParamsNH.mat'); 
    load('estimates/ParamsH.mat'); 

    ParamsH.beta = ParamsH.alpha ;
    
    % indicator for value inclusive of tax/subsidy
    ParamsNH.tax = 1 ;
    ParamsH.tax  = 1 ;

    %% 1/7 baseline decentralized equilibrium

    disp(' ')
    disp('+++++++++++++++++++++++++++++++++')
    disp('+++ 1/7 Computing baseline DE +++')
    disp('+++++++++++++++++++++++++++++++++')
    disp(' ')

    options.Hosios = 0 ;
    options.IdenticalLocations = 0 ;
    options.Policy = 'Beta' ; % 'Beta' or 'Alpha'

    Params = ParamsNH ;
    run AllModelParameters.m
    
    % for possible policy exercise
    Grids.DlogPolicy = zeros(Grids.Nx,1) ;
    Grids.tax        = ones(Grids.Nx,1) ;

    Params.ineff = 1 ;
    run BargainingPower.m

    GE.vU     = 0.0626 ;     % initial guess
    GE.vL     = 0.0036 ;     % initial guess
    GE.C      = 9.3403 ;     % initial guess
    Params.ce = 5.6072e-06 ; % initial guess

    [ vC, zetaC , LhC , LC , vUC , CC , ceC , ...
      popC , fC , wC , whC , uC , HC , hC , sC ] = SolveEntry( Grids , Params , GE ) ;

    GEC.vL = vC(1) ;
    GEC.vU = vUC ;
    GEC.C = CC ; 

    ParamsC = Params ;
    GridsC  = Grids ;
    ParamsC.ce = ceC ;

    run ComparisonDE.m       % outputs objects for ZFU exercise


    %% 2/7 planner case with decentralized equilibrium estimation

    disp(' ')
    disp('++++++++++++++++++++++++++++++++++++++++++++++')
    disp('+++ 2/7 Computing planner w/ DE parameters +++')
    disp('++++++++++++++++++++++++++++++++++++++++++++++')
    disp(' ')

    options.Hosios = 1 ;
    options.IdenticalLocations = 0 ;
    options.Policy = 'Alpha' ; % 'Beta' or 'Alpha'

    Params = ParamsNH ;
    run AllModelParameters.m
    
    % for possible policy exercise
    Grids.DlogPolicy = zeros(Grids.Nx,1) ;
    Grids.tax        = ones(Grids.Nx,1) ;

    Params.ineff = 0 ;
    run BargainingPower.m

    Params.ce = ceC ; % actual value
    
    GE.vU = 0.0572 ; % initial guess
    GE.vL = 0.0280 ; % initial guess
    GE.C  = 6.4221 ; % initial guess

    [ vP, zetaP , LhP , LP , vUP , CP , ceP , ...
      popP , fP , wP , whP , uP , HP , hP , sP ] = ...
        SolveEq( Grids , Params , GE ) ;

    GEP.vU = vUP ;
    GEP.vL = vP(1) ;
    GEP.C = CP ;

    ParamsP = Params ;
    GridsP  = Grids ;

    ParamsP.ce = ceP ;

    % labor subsidy that would implement the planner's allocation
    run LaborSubsisdy.m


    %% 3/7 planner case with planner estimation

    disp(' ')
    disp('+++++++++++++++++++++++++++++++++++++++++++++++++++')
    disp('+++ 3/7 Computing planner w/ planner parameters +++')
    disp('+++++++++++++++++++++++++++++++++++++++++++++++++++')
    disp(' ')

    options.Hosios = 1 ;
    options.IdenticalLocations = 0 ;
    options.Policy = 'Alpha' ; % 'Beta' or 'Alpha'

    Params = ParamsH ;
    run AllModelParameters.m
    
    % for possible policy exercise
    Grids.DlogPolicy = zeros(Grids.Nx,1) ;
    Grids.tax        = ones(Grids.Nx,1) ;

    Params.ineff = 0 ;
    run BargainingPower.m


    Params.ce = 4.9764e-07 ; % initial guess
    GE.C      = 476.1297 ;   % initial guess

    [ vPP, zetaPP , LhPP , LPP , vUPP , CPP , cePP , ...
      popPP , fPP , wPP , whPP , uPP , HPP , hPP , sPP ] = SolveEntry( Grids , Params , GE ) ;

    GEPP.vU = vUPP ;
    GEPP.vL = vPP(1) ;
    GEPP.C = CPP ;

    ParamsPP = Params ;
    GridsPP  = Grids ;

    ParamsPP.ce = cePP ;


    %% 4/7 quasi-optimal policy (decentralized equilibrium w/ policy)

    disp(' ')
    disp('++++++++++++++++++++++++++++++++')
    disp('+++ 4/7 Quasi-optimal policy +++')
    disp('++++++++++++++++++++++++++++++++')
    disp(' ')

    options.Hosios = 0 ;
    options.IdenticalLocations = 0 ;
    options.Policy = 'Beta' ; % 'Beta' or 'Alpha'

    Params = ParamsNH ;
    run AllModelParameters.m
    
    % for possible policy exercise
    Grids.DlogPolicy = zeros(Grids.Nx,1) ;
    Grids.tax        = tauVA ; % tauVAadjInt ;

    Params.ineff = 1 ; % combined with Policy = Beta, get bargaining power = beta, 
                       % as if correct the inefficiency in the location
                       % decision of firms
                       
    Grids.tax        = tauVA ;
    Grids.DlogPolicy = [ log(tauVA(end))   - log(tauVA(end-1)) ; ...
                         log(tauVA(2:end)) - log(tauVA(1:end-1)) ] ./ Grids.dx ; 

    run BargainingPower.m

    GE.vU = 0.0585 ;         % initial guess
    GE.vL = 0.0397 ;         % initial guess
    GE.C  = 10.2182 ;        % initial guess
    Params.ce = ParamsC.ce ; % actual value

    [ vQO, zetaQO , LhQO , LQO , vUQO , CQO , ceQO , ...
      popQO , fQO , wQO , whQO , uQO , HQO , hQO , sQO ] = SolveEq( Grids , Params , GE ) ;

    GEQO.vL = vQO(1) ;
    GEQO.vU = vUQO ;
    GEQO.C = CQO ; 

    ParamsQO = Params ;
    GridsQO  = Grids ;


    %% 5/7 identical locations

    disp(' ')
    disp('+++++++++++++++++++++++++++++++')
    disp('+++ 5/7 Identical locations +++')
    disp('+++++++++++++++++++++++++++++++')
    disp(' ')

    options.Hosios = 0 ;
    options.IdenticalLocations = 1 ;
    options.Policy = 'Beta' ; % 'Beta' or 'Alpha'

    Params = ParamsNH ;
    run AllModelParameters.m
    
    % for possible policy exercise
    Grids.DlogPolicy = zeros(Grids.Nx,1) ;
    Grids.tax        = ones(Grids.Nx,1) ;

    Params.ineff = 1 ;
    run BargainingPower.m

    GE.vU = 0.0452 ;  % initial guess
    GE.vL = 0.0090 ;  % intial guess
    GE.C  = 7.4683 ;  % initial guess
    Params.ce = ceC ; % actual value

    [ vI, zetaI , LhI , LI , vUI , CI , ceI , ...
      popI , fI , wI , whI , uI , HI , hI , sI ] = SolveEq( Grids , Params , GE ) ;

    GEI.vL = vI(1) ;
    GEI.vU = vUI ;
    GEI.C = CI ; 

    ParamsI = Params ;
    GridsI  = Grids ;

    ParamsI.ce = ceI ;


    %% 6/7 ZFU

    disp(' ')
    disp('++++++++++++++++')
    disp('+++ 6/7 ZFUs +++')
    disp('++++++++++++++++')
    disp(' ')

    options.Hosios = 0 ;
    options.IdenticalLocations = 0 ;
    options.Policy = 'Beta' ; % 'Beta' or 'Alpha'

    Params = ParamsNH ;
    run AllModelParameters.m

    run PrepareZFU.m
    % Measure of locations in ZFU = 0.064408
    % ZFU subsidy = 19.4331%

    Params.ineff = 1 ;
    run BargainingPower.m

    GE.vU = GEC.vU  ; % initial guess
    GE.vL = GEC.vL  ; % initial guess
    GE.C  = GEC.C   ; % initial guess
    Params.ce = ceC ; % actual value

    [ vZ, zetaZ , LhZ , LZ , vUZ , CZ , ceZ , ...
      popZ , fZ , wZ , whZ , uZ , HZ , hZ , sZ ] = SolveEq( Grids , Params , GE ) ;

    GEZ.vL = vZ(1) ;
    GEZ.vU = vUZ ;
    GEZ.C = CZ ; 

    ParamsZ = Params ;
    GridsZ  = Grids ;
    
    run Redistribution.m
    

    %% 7/7 alternative parameters

    disp(' ')
    disp('++++++++++++++++++++++++++++++++++')
    disp('+++ 7/7 alternative parameters +++')
    disp('++++++++++++++++++++++++++++++++++')
    disp(' ')

    options.Hosios = 0 ;
    options.IdenticalLocations = 0 ;
    options.Policy = 'Beta' ; % 'Beta' or 'Alpha'

    Params = ParamsNH ;
    run AllModelParameters.m
    
    % for possible policy exercise
    Grids.DlogPolicy = zeros(Grids.Nx,1) ;
    Grids.tax        = ones(Grids.Nx,1) ;

    Params.ineff = 1 ;
    run BargainingPower.m

    % assign new parameters
    Params.gamma = 3 ; 

    GE.vL = 2.8639e-06 ; % initial guess
    GE.vU = 3.5490e-04 ; % initial guess
    GE.C  = 0.0108 ;     % initial guess
    Params.ce = ceC ;    % actual value

    [ vA, zetaA , LhA , LA , vUA , CA , ceA , ...
      popA , fA , wA , whA , uA , HA , hA , sA ] = SolveEq( Grids , Params , GE ) ;

    GEA.vU = vUA ; 
    GEA.vL = vA(1) ; 
    GEA.C  = CA ; 
    
    ParamsA = Params ;
    GridsA  = Grids ;

    if options.Save
        save('estimates/JMP.mat')
    end


end



%% results

% load outcomes
if options.Load
    load('estimates/JMP.mat') ;
end

% saving results (graphs, tables)
options.SaveResults = 0 ;

% compute all moments
MomentsCE = ComputeMoments( ParamsC  , GridsC  , LC  , LhC  , zetaC  , vC  , GEC  ) ;
MomentsQO = ComputeMoments( ParamsQO , GridsQO , LQO , LhQO , zetaQO , vQO , GEQO ) ;
MomentsPP = ComputeMoments( ParamsPP , GridsPP , LPP , LhPP , zetaPP , vPP , GEPP ) ;
MomentsZ  = ComputeMoments( ParamsZ  , GridsZ  , LZ  , LhZ  , zetaZ  , vZ  , GEZ  ) ;
MomentsI  = ComputeMoments( ParamsI  , GridsI  , LI  , LhI  , zetaI  , vI  , GEI  ) ;

% load moments in data
run LoadDataDecomp.m

% plot of JL and JF against u/(1-u) in DE
run FigureScatter.m

% table: different estimations and data
run ComparisonEstimations.m

% table: compare both different estimations, no pooling (QO), identical locations and data
run ComparisonDEversions.m

% figure: allocations under DE and policies
run FigurePolicies.m

% table: welfare gains from policies
run AgWelfareGains.m

% figures: welfare gains from policies
run WelfareGainsLocation.m

% regressions of JL, JF and U on wages and population
run WagePopReg.m

% ZFU DiD regressions
run DiDZFU.m





